library(mvtnorm)
library(tmvtnorm)
library(bayesm)
dyn.load("test/Orthant.so")

source("R/test_Orthan_S.R")
Sub<-function(A, i ,j)
{
	p<-nrow(as.matrix(A))	
	if(i>j)
	{
		s<-j
		j<-i
		i<-s
	}
	Y<-matrix(.C("Substitute",as.double(A),res=as.double(A),as.integer(p),as.integer(i),as.integer(j))$res,nrow=p)
}
order<-function(S,a,m,M)
{
	pp<-ncol(as.matrix(S))
	a<-a-m
	for(i in 1:pp)
	{	
		resp<-1
		for( j in i:pp)
		{
			SS<-Sub(S,i,j)
			SS<-SS[1:i,1:i]
			aa<-a
			aa[i]<-a[j]
			aa<-aa[1:i]
			p<-ncol(as.matrix(SS))
			L<-t(chol(SS))
			theta<-0
			res<-.C("ghk_order", as.double(L),as.double(aa),as.double(rep(0,p)),as.integer(p),res=as.double(theta))$res
			if(res<resp)
			{
				resp<-res
				k<-j
			}		
		}
			S<-Sub(S,i,k)
			s<-a[i]
			a[i]<-a[j]
			a[j]<-s
	}
	return(list(S=S,a=a))
}

jt<-function(S,a,m,M,nu)
{
	p<-ncol(S)
	L<-t(chol(S))
	theta<-0
	tr<-1
	res1<-.C("sampler5", as.integer(M),as.integer(p), as.double(L),as.double(a-m), as.integer(nu),as.double(0.9), res=as.double(theta))
	return(res1$res)
}

mt<-function(S,a,m,M,nu)
{
	p<-ncol(S)
	L<-t(chol(S))
	theta<-0
	res1<-.C("sampler6", as.integer(M),as.integer(p), as.double(L),as.double(a-m),as.integer(nu),as.double(rep(0,p)),as.double(0.9), as.integer(1),res=as.double(theta))
	return(res1$res)
}


library(MASS)
data(Pima.tr)
library(arules)

X<-as.matrix(Pima.tr[1:7])
source("~/Documents/Memoire/src/Generate.R")
source("~/Documents/Memoire/src/new_data.R")
X<-X[,-1]
M<-500
p<-10
nu<-10
#set.seed(3)
A<-matrix(rcauchy(p*p,0.001),nrow=p)
#source("R/Covariance.R")
A<-t(A)%*%A
 
g<-rep(0,3)
a<-rep(0.1,p)
m<-rep(0,p)

t<-proc.time()
reso<-order(A,a,m,M)
tOrd<-(proc.time()-t)

cat("i")
t<-proc.time()
resj<-jt(A,a,m,M,nu)
h<-(proc.time()-t)

t<-proc.time()
resj<-jst(A,a,m,M,nu)
jsti<-(proc.time()-t)
		
cat("i")
t<-proc.time()
b<-jt(reso$S,reso$a,m,M,nu)
jj<-(proc.time()-t)+tOrd

cat("i")
t<-proc.time()
solm<-mt(reso$S,reso$a,m,M,nu)
mi<-proc.time()-t+tOrd

cat("i")
t<-proc.time()
solm<-mt(A,a,m,M,nu)
j<-proc.time()-t

Mh<-floor(j[[3]]/h[[3]]*M)
Mjj<-floor(j[[3]]/jj[[3]]*M)
Mmi<-floor(j[[3]]/mi[[3]]*M)
Mjst<-floor(j[[3]]/jsti[[3]]*M)

cat("\n")
cat("\n")
cat(resj)
cat("\n")
cat(res3)
#source("R/test_Orthan_S.R")
#solst<-jst(A,a,m,M/10,1,1,10)
#
Mm<-10
#XX<-rtmvnorm(10000,mean=m,sigma=A,lower=a,upper=rep(Inf,p))
i=1
b<-matrix(0,nrow=Mm,ncol=5)
for(i in 1:Mm)
{
	cat("i")
	#b[i,1]<-jt(A,a,m,Mh,nu)
	b[i,2]<-jt(reso$S,reso$a,m,Mjj,nu)
	b[i,3]<-mt(reso$S,reso$a,m,Mmi,nu)
	b[i,4]<-mt(A,a,m,M,nu)
	b[i,5]<-jst(reso$S,reso$a,m,M,nu)
	
	#write.table(b,"res_Orthan.txt")
	cat("loop")
	cat(i)
	cat("\n")
}
#.C("testbrt")
#pdf("ghkgibbs_RD.pdf")
z<-factor(c(rep("sep",Mm),rep("I-Res",Mm),rep("Res",Mm),rep("Mix",Mm)))
#z<-factor(c(rep("No-Res",Mm),rep("I-Res",Mm),rep("I-NoRes",Mm),rep("I-Res_RD",Mm)))
boxplot(as.vector(exp(b[,2:5]))~z)
#dev.off()
dyn.unload("test/Orthant.so")
#dyn.unload("test/RLogit.so")
#YY<-rmvnorm(10000,mean=apply(XX,2,mean),sigma=var(XX))
#plot(YY)
#points(XX,col="red")
#dyn.unload("test/Rprobit.so")
#dyn.unload("test/IBIS.so")
